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Abstract. 

We present a new numerical method for the generation of binary neutron star initial 
data using a method along the lines of the the Wilson-Mathews or the closely related 
conformal thin sandwich approach. Our method uses six different computational 
domains, which include spatial infinity. Each domain has its own coordinates which 
are chosen such that the star surfaces always coincide with domain boundaries. These 
properties facilitate the imposition of boundary conditions. Since all our fields are 
smooth inside each domain, we are able to use an efficient pseudospectral method 
to solve the elliptic equations associated with the conformal thin sandwich approach. 
Currently we have implemented corotating configurations with arbitrary mass ratios, 
but an extension to arbitrary spins is possible. The main purpose of this paper is to 
introduce our new method and to test our code for several different configurations. 
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1. Introduction 

Currently several gravitational wave detectors such as LIGO 2j, Virgo ^ ^4] or 
GEO [5] are already operating, while several others are in the planning or construction 
phase [6]. One of the most promising sources for these detectors are the inspirals and 
mergers of binary neutron stars (NS). In order to make predictions about the final phase 
of such inspirals and mergers, fully non-linear numerical simulations of the Einstein 
Equations are required. To start such simulations initial data are needed. The emission 
of gravitational waves tends to circularize the orbits [T, '8]. Thus, during the inspiral, 
we expect the two NSs to be in quasicircular orbits around each other with a radius 
which shrinks on a timescale much larger than the orbital timescale. This means that 
the initial data should have an approximate helical Killing vector In addition, one 
would like to have the initial data in coordinates such that this approximate symmetry is 
manifest, i.e. the time evolution vector should lie along so that the time derivatives 
of the evolved quantities are minimized. In order to achieve these goals we use the 
Wilson-Mathews approach [9l[T0], which is closely related to the conformal thin sandwich 
formalism [Tl]. The Wilson-Mathews approach approach has already been successfully 
used by several groups. Among them are results for corotating [121 lEl [IH [E] and 
irrotational |16l [T71 [181 EH [20] NS binaries with equal masses. One group has also 
produced results for unequal mass systems [21], [22] . 

In this paper we present a new numerical method to construct initial data for binary 
NSs in corotating configurations for arbitrary mass ratios. The main focus is on the 
numerical method rather than on new physics. We describe an efficient implementation 
of this method with the SGRID code [23], discuss code tests, and compare with previous 
results. 

Throughout we will use units where G = c = 1. Later when we present numerical 
results we will use fully dimensionless units by setting k in the polytropic equation of 
state to K = G = c = 1. Latin indices such as i run from 1 to 3, while Greek indices 
such as /i run from to 3. The paper is organized as follows. In Sec. [2] we describe the 
General Relativistic equations that govern binary neutron stars described by perfect 
fluids. Sec. [3] describes our particular numerical implementation of these equations, 
followed by results for some particular configurations in Sec. [H We conclude with a 
discussion of our method in Sec. [5] 

2. Binary neutron stars in General Relativity 

In this section we describe the equations governing binary NSs in quasi circular orbits. 
2.1. ADM decomposition of Einstein's equations 

We use the Arnowitt-Deser-Misner (ADM) decomposition of Einstein's equations (see 
e.g. [21]) and write the line element 

ds^ = -a^dt^ + jij{dx' + (3'){dx^ + (3^) (1) 
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in terms of the lapse a, shift /3* and the 3- metric 7jj. The extrinsic curvature is defined 
by 

Kij = -^{dtlij - £plij), (2) 
Za 

With these definitions Einstein's equations spht into the evolution equations 

dtlij = - 2aKij + ffs'jij 

dtKi, = aiRij - 2KuK\ + KK^^) - D.D^a + f^K^j 

-8nS^j+A7c^ij{S-p) (3) 
and the Hamiltonian and momentum constraint equations 
R - KijK'^ + K'^ = IQnp 

Dj{K'^-^'^K) =87Tj\ (4) 

Here Rij and R are the Ricci tensor and scalar computed from ■jij, Di is the derivative 
operator compatible with and all indices here are raised and lowered with the 3-metric 
7jj. The source terms p, i\ Sij and S = Sij are projections of the stress-energy tensor 
T^^ given by 

f = - T^^n^r 

S'' = T^u^Y' (5) 

and correspond to the energy density, flux and stress-tensor. The vector appearing 
here is the the 4-vector normal to a t = const slice. 

2.2. Decomposition of 3-metric and extrinsic curvature 

As in [9l [To] the 3-metric 7,^ is decomposed into a conformal factor ip and a conformal 
metric 7,^ such that 

7ij = i^%3- (6) 

The extrinsic curvature is split into its trace K and its tracefree part Aij by writing it 

as 

K^J = A, + ^7,,ir (7) 

2.3. Quasi equilibrium assumptions 

We now make some additional simplifying assumptions. First we assume that our binary 
is in an approximately circular orbit and that the stars are corotating. This implies the 
existence of an approximate helical Killing vector In a coordinate system where 
this helical symmetry is manifest and the time evolution vector lies along all time 
derivatives should approximately be zero. Here we only assume that the time derivative 
dt'^ij of the conformal metric and the time derivative dtK of the trace of the extrinsic 
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curvature vanish. The former allows us to express the extrinsic curvature in terms of 
the shift and results in 

where 

(L/Sy^ = &I3^ + &(3' - '^D,f3\ (9) 

and Dk is the derivative operator compatible with 7jj. The assumption dtK = together 
with the evolution equation of K (derived from Eq. ([3])) implies 

ij-^[DkD''{axlj) - aDkD^tl)] = a{R + Kf + p'DiK 

+ 47ra(^-3p). (10) 

2.4- Further simplifications and boundary conditions 

Next we also choose a maximal slice and thus K = and assume that the conformal 
3-metric is flat and given by [9|, [10] 

lij = Sij. (11) 

This latter assumption merely simplifies our equations and could in principle be 
improved by e.g. choosing a post-Newtonian expression for as in [251 [26]. Using 
Eq. (JTTil the Hamiltonian and momentum constraints in Eq. and Eq. (ITOl) simplify 
and we obtain 

D'^P = - ^{LBriLBy, - 2n^'p 

D\a^) =a^ -^{LBy^{LBy,+2TT^p\p + 2S) 

(12) 

where {LbY^ = D'B^ + B' - p'WkB'', Di = di, and 

B' = (3' + uje'^%x^ -x'cm)- (13) 

Here x'^^j^ denotes the center of mass position and u is the orbital angular velocity, 
which we have chosen to lie along the ^-direction. The elliptic equations ([T2!) have to 
be solved subject to the boundary conditions 

lim ^ = 1, lim B' = 0, lim a^b = 1. (14) 

r— +00 r—>oo r—>oo 

at spatial infinity. 
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2.5. Matter equations 

We assume that the matter in both stars is a perfect fluid with a stress-energy tensor 
T^- = [po(l + e) + P]m^m- + Pg>^r (15) 

Here Pq is the mass density (which is proportional the number density of baryons) , P is 
the pressure, e is the internal energy density divided by po, is the 4-velocity of the 
fluid and g^'^ is the spacetime metric. The matter variables in Eq.([5]) are then 

p =a2[po(l + e)+P]MV-P 

= «[po(l + e) + P]mV(m7m° + 13') 
S'^ = [po(l + e) + P]mV(m7u° + (3'){u^/u^ + (3^) 

+ PY^ (16) 
The fact that V^T^^ = yields the relativistic Euler equation 

[po(l + e) + PW.u^ = -{g^'' + u^u'')V.P, (17) 
which together with the continuity equation 

^u{poun = (18) 
governs the fluid. 

For corotating stars we can show that the continuity equation ( fT8|) is identically 
satisfied. Furthermore one can show that the Euler equation leads to (see e.g. problem 
16.17 in [27]) 

[p^{l + e) + P]d\n{u^e) = -dP. (19) 
where is the assumed helical Killing vector. With the help of the first law of 
thermodynamics ((i[po(l + e)] = [po(l + e) + P]dpo/po) this equation can be integrated 
to yield 

"po(l + e) + P' ^'"^ 
where C1/2 are constants of integration for each star. We will later choose them such 
that the rest mass of each star has a prescribed value. In corotating coordinates and 
taking into account our conformally fiat 3-metric, m^^^ can be written as 

= -l/n° = -[a' - (21) 

In order simplify the problem we assume a polytropic equation of state 

P = KpJ+'/". (22) 

It is then convenient to introduce the dimensionless ratio 

q = Pjp^, (23) 

which we use to write 

Po = 
P = 

e = nq. (24) 
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3. Numerical method 



In order to construct binary NS initial data we have to solve the five elliptic equations 
in Eq. f|T2l) . with the matter terms given by Eqs. f|T6|) . fl2T|) and fl24l) . In addition, our 
data also have to satisfy Eq. ( !20l) . which can be expressed as 

1 / C'l/2 



n + 1 \u^ii^ 



1 



(25) 



for each star. We will solve the whole set of equations by iterating over the following 
steps: (i) We first come up with an initial guess for q in each star, in practice we simply 
choose Tolman-Oppenheimer-Volkoff (TOV) solutions (see e.g. Chap. 23 in [21]) for 
each, (ii) Next we solve the 5 coupled elliptic equations (fT2l) for this given q. (iii) Then 
we use Eq. (!25|) to update q in each star. The constants C1/2 in general have different 
values for each star. We adjust the value for each star such that it has a prescribed 
rest mass. After updating q we go back to step (ii) and iterate until all equations are 
satisfied up to a given tolerance. 



3.1. Coordinates adapted to star surfaces 

Note that the matter is smooth inside the stars. However, at the surface (at g = 0), po? 
P and e are not differentiable. This means that if we want to take advantage of a spectral 
method, the star surfaces should be domain boundaries. A difficulty with our iterative 
approach, however, is that each time we update q the matter distributions change, so 
that the stars change shape or even move. Hence the domain boundaries have to be 
changed as well. In order to address this problem we introduce several domains each 
with its own coordinates. These coordinates depend on two freely specifiable functions 
which will allow us to vary the location of the domain boundaries, so that we can always 
adapt our domains to the current star surfaces in each iteration. As in the initial data 
approaches in [28l [29] our aim was to introduce as few domains as possible. 

The coordinates we will use, are very similar to the ones introduced by 
Ansorg [30]. We place both stars on the x-axis and write down the necessary coordinate 
transformations in two steps. First we express the standard Cartesian coordinates as 



X 



(X2 + /?2)2 

1 

(X2 + i?2)2 
1 



+ 1 



(X^ - 



XRcos i 



(X2 + i?2)2 ^ ^^^i^'^' (26) 

where 6 is a parameter related to the distance between the stars, and X, R are functions 
of the new coordinates [A, B, (p) we will use in each domain. Note that spatial infinity 
is located at the point where X = R = and that in order to cover all {x, y, z) it 
is sufficient the restrict X,R,(j) to the ranges 0<X<l,0<i?< yl — X^ and 
< < 27r. In order to complete the coordinate transformation between the Cartesian 
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{x, y, z) and the new coordinates (A, 5, 0), we now write down X and R as functions of 
(A, 5,0). Inside starl we use 

X = (1 - A){3?[C+(i?, 0)] - 0)]} 

+ 5 cos([l - A] arg[C+(l, 0)]) + (1 - B)A 
R = {1- A){^[C+{B, 0)] - B^[C^{1, 0)]} 

+ 5sin([l-A]arg[C+(l,0)]), (27) 

where the strictly positive function o"+(i?,0) in 



\ 



f MM)±M) (28) 



determines the shape of the star surface. The surface is always located at A = 
but depending on the choice for o"+(i?,0) it will be at different {x,y,z), e.g. for 
(j+(i?,0) = const we will get a spherical surface in {x,y,z). Note that this star is 
located around x = b. Inside star2 we use a similar transformation given by 

X = (1 - A){^[C4B, 0)] - B^[C41, 0)]} 
+ Bcos{^A+[l-A] arg[C_(l,0)]) 



R = (l-A){53[C_(i?,0)]-i?$5[C_(l,0)]} 

+ B sin(^A +[1-A] arg[C_(l, 0)]) + (1 - B)A, 



(29) 



t.nJ'-^^4±^] (30) 



where the strictly negative function o'-{B,(j)) in 

C_(5,0) = ^ ^ ^ J 

determines where the star surface {A = 0) is located in {x, y, z) coordinates. Star2 is 
located around x = —b. Note that the A,B,(j) coordinates are different inside each star, 
but in order to cover each star their ranges are 0<y4<l,0<i?<l and < < 27r 
in each case. 

The outside of both stars is covered by two additional domains. The first one covers 
the region outside starl for all positive x, while the second one covers the region outside 
star2 for all negative x. Both coordinate transformations can be written as 

X = (1 - A){^[C±iB, 0)] - 53?[C±(1, 0)]} 

+ Bcos{^A+[l-A] arg[C±(l,0)]) 
R = {1- A){^[C±{B, 0)] - B^[C^{1, 0)]} 

+ B sin(^A +[1-A] arg[C±(l, 0)]), (31) 

where we use C+{B, 0) in the former and C-{B, 0) in the latter. In each case the star 
surface is at A = and spatial infinity is at {A, B) = (1, 0). 
Figure [1] shows the coordinate lines in z = plane. 
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(A,B)=(0,0) 




A-0 



(A,B)-(0,1) 



Figure 1. The plot shows the hnes of constant A and B in the xy-plane for b — 1.84, 
a+{B,<j>) = 1.51 and a-{B,(f) = —1.28. The two neutron stars are marked with NSl 
and NS2. In addition, the (a;, y) positions of a few points are indicated. The coordinate 
singularities at A = 1 inside the stars are located ai x — ±.b. 



3.2. Spectral method 



In order to solve the elliptic equations f|T2l) we use the SGRID code |23j which employs 
pseudospectral methods to accurately compute spatial derivatives. We use Chebyshev 
expansions in the A- and 5-directions and Fourier expansions 0-direction. As collocation 
points we choose 

/ TXI 



A, 



1 
2 
1 
2 



cos 



1 — cos 




(32) 



where /, j, k are integers obeying 

0<l KriA, 0<j<nB, 0<k<n^. (33) 

The number n^, and of collocation points in each direction is chosen to be equal 
in all four domains, to ensure that the grid points on the boundaries of two adjacent 
domains will be at the same {x, y, z) location. As in [23] will solve Eq. (fT2|) as written 
down in Cartesian form and compute derivatives like dxil' using the chain rule: 



ox ox OX 



(34) 



Note that all points with with B = Q oi B = 1 lie along the x-axis, with x 
independent of 0. Hence along this axis we have the standard coordinate singularity 
of polar coordinates. Furthermore, all points with A = 1 in the interior of each star 
correspond to just one point on the x-axis. Thus there is an additional coordinate 
singularity aX A = 1 inside each star. We have found that if we simply use the A, B, (j) 
coordinates as described above we were not able to solve the elliptic equations (fT2|) . 
The culprit is the singularity ai A = 1 inside each star. Near these points the Jacobian 
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Eq. flM|l . This problem would also occur if we did use the exact same coordinates as 
proposed by Ansorg [30]. In [30] the problem is not addressed since for black holes one 
can use excision boundary conditions and does not need the inner domains. One way 
around this problem would be to construct different basis functions, which would have 
to be chosen such that they have vanishing B and derivatives ai A = 1. However, then 
we would not be able to use Fast Fourier transforms anymore to compute derivatives. 
For this reason we have chosen a different approach. We simply restrict the range of A 
inside each star so that inside < A < A^ax < 1- The collocation points in A inside 
the stars are then given by 



We typically choose Amax = 0.85. In this way we completely avoid the singularities at 
A = 1. Of course then, our inner domains do not cover the entire star interiors any 
more. Instead they leave out a small hole around A = 1. We simply cover this hole by 
placing two additional cubical domains inside each star. Each cube is chosen such that 
completely covers the hole (described by A^ax < ^ < !)• In each cube we use standard 
Cartesian coordinates so that it overlaps with part of the inner domain covered by the 
A, B, (J) coordinates. The collocation points inside each cube are then 



where the minima and maxima Xmin, x^ax, Vmin, Umax Zmin, z^ax are chosen such that 
the cube covers the hole, and where nc = 6 is typically sufficient, since the cubes are 
very small. Figure [2] shows how a cube is fitted into starl. 

In order to numerically solve the elliptic equations in Eq. ( |T2l) we arrange the values 
of the fields B^ and a%l) at each grid point in a vector w. Since we solve for 5 fields, 
the dimension of this vector is five times the total number of grid points. In order to find 
a w that satisfies Eq. (fT2l) we impose Eq. ( JT2l) at all interior grid points. At adjacent 
domain boundaries we impose the conditions that fields and their normal derivatives 
are equal on both sides. At infinity we impose Eq. ([HI). In the domains covered by the 
A, 5, coordinates we impose the following regularity conditions along the x-axis: For 
> we demand 




(35) 




(36) 



^{AuBj.cPk) = ^(Aj,S„0o) 



(37) 



while for /c = we impose 



ds^{Ai, Bj, 0o) + dsd^d^^iAi, Bj, 0o) = 0. 



(38) 
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Figure 2. In order to avoid the coordinate singularity at A = 1 inside the stars, we 
have restricted the range of ^ to < A < 0.85. The remaining space is fiUed with 
a small cube of length 0.0166. This figure is a blowup of the inside of NSl in Fig. [1] 
where the small cubes are not visible. 

Here \E' stands for either ip, or aip, and s = ^ is the distance from the x- 
axis. In order to deal with the cubes which overlap the A, 5, covered domains inside 
each star, we impose the condition that the fields at points on the cube boundaries 
must be equal to the fields in the A,B,(f) covered domain interpolated to these points. 
Correspondingly we also demand that the fields at points on the A = Amax boundaries 
must equal to the fields in the cube interpolated to these points. This interpolation 
is done with our given spectral accuracy. To compute a field at a point that is not 
a collocation point, we first compute the spectral expansion coefficients from the field 
values at the collocation points in the domain of interest. To interpolate to any point 
in this domain we then compute the field value from a sum over coefficients times the 
basis functions evaluated at the point in question. This last interpolation step can 
be computationally expensive if we interpolate onto many points, while going between 
field values and collocation points can be done via Fast Fourier transforms and is thus 
not very expensive. Note, however, that in our case these interpolations are not too 
costly, because our cubical domains have only 6^ grid points. Such a small number of 
points should always be sufficient because the cubic domains are so small so that all 
fields are nearly constant inside the cubes, e.g. on the scale of Fig. 1 the cubes are 
not visible. This is an advantage over the domain decomposition used in [311 [T7j where 
interpolations are needed between domains with many more grid points. 

If we take all these conditions into account we obtain N = 'inAnBn^ + 2n^^ non- linear 
equations of the form 

/^H = 0, m = l,2,...,iV (39) 

for the N unknowns comprising the solution vector w. We solve this system of equations 
by a Newton- Raphson scheme. This scheme requires an initial guess, for which we simply 
use two TOV solutions in conformally flat isotropic coordinates. In order to solve the 



A new numerical method to construct binary neutron star initial data 



11 



linearized equations 

= -/.»(-) («) 

in each Newton-Raphson step, we note that fm{w) contains spectral derivatives of w 
in different directions, so that the N x N matrix ^'^iT'' is sparse in the sense that it 
contains about 95% zeros. So in order to numerically solve the linearized Eq. (HOj) we 
use the sparse matrix solver UMFPACK 



3.3. Iteration scheme 

As already mentioned solving the elliptic equations in Eq. flT2|) once is not enough. After 
each solve, we have to adjust q using Eq. ( l25l) . This adjustment presents the problem 
that the star surfaces (located at g = 0) change. Hence we have to also adjust 0) 
and o"_(i?, (f)) in order to keep both star surfaces at A = 0. While this adjustment is not 
hard to implement, it incurs a high computational cost since adjusting a+{B,(f)) and 
a_{B,(j)) amounts to changing our computational grid, and after each such adjustment 
we need to interpolate all relevant fields onto the new grid. In addition, the adjustment 
has to be carried out several times after each individual elliptic solve, since we were only 
able to achieve a stable iteration scheme if we pick the free constants Ci, C2, OJ and 
xcM as follows. Let us call the intersection points of the x-axis with the side of the star 
surface not facing the origin Xouti and Xout2 (located at {A,B) = (0,0)) for each star. 
We then determine u and xcm by requiring that Xouti and Xout2 remain constant. This 
task is accomplished by a root finder. In each iteration of this root finder Ci, C2 are 
adjusted such that the rest mass of each star remains constant (by another root finder). 
These root finders have to evaluate q and thus need to adjust the domain shapes several 
times. If we do not adjust u and xcm we find that the stars drift around too much for 
the iterations to converge. Instead of Xouti and Xout2 it is also possible to fix the points 
Xmaxi and Xmax2 wherc the maximum values of q occur. 

In addition, we have observed that once a new q is set the elliptic solve often 
"over corrects" which again can result in an unstable iteration scheme. To overcome this 
problem we typically do not take the ip, and aip coming from solving Eq. ( |T2l) as our 
new fields. Rather, we take the average of this solution and the and aip from the 

previous iteration step as our new fields. In this way ip, B^ and aip change less from 
one iteration step to the next. 



4. Results 

All the results presented in this section were computed for n = 1 polytopes (see Eq. (l22l) ) 
in units where k = 1. 

In order to check that our SGRID code is working properly we have checked the 
convergence of the constraints. For these tests we have computed the constraints 
directly from Eq. (jl]) for different numbers of grid points. From Fig. [3] we see that 
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Figure 3. The L^-norm of the Hamihoniaii constraint on the grid inside star 1 
converges exponentiaUy with the number of grid points ^ hb = n. The plot is for 
an equal mass binary with rest masses moi = mo2 = 0.05952. The star centers are 
located at Xmaxi = —Xmax2 — 2.112, and = 8 and rtc = 6 are kept fixed. 

the Hamiltonian constraint converges exponentially with the number of grid points, as 
expected for a spectral method. The momentum constraints as well as dtK converge to 
zero in a similar fashion. 

We have computed initial data for the configurations listed in table [H Each 
configuration is described by the rest masses mgi and of the two stars given by 

moi=/ powW^^'a;, ^ = {1,2} (41) 
Jstari 

and the separation parameter 6, which appears in the coordinate transformations and 
is approximately half the separation. For each configuration we have also computed the 
ADM mass and angular momentum given by 

Madm = j {p + ■^^2(^Bf(^B)ii) ^'^'^ (42) 

and 

Jadm = Jli^- ^CM)f - yf] i^'^d^x. (43) 

Our iterative scheme also yields the orbital angular velocity uj and the location of the 
center of mass xcm- Furthermore we list the maximum values of q in each star along the 
X-axis, together with their x-coordinates, and also the locations of the inner and outer 
edges of each star. We also show the distance di2 between the stars and star diameters 
di/2 defined by di2 = \xmaxi - Xmax2\ and di/2 = \xouti/2 - a;i„i/2|- The equal mass 
configuration in table [U is very close to configurations already computed by Baumgarte 
et al. [13] and also by Gourgoulhon et al. [T7] and agrees with those to better than 
1% (see Table II. in [H]). Note that our code has no problems handling unequal mass 
systems, as well as systems that are far apart. In addition, it is very memory efficient. 
A typical run with = = 18 points needs only about 80MB of memory. The initial 
data can thus be generated on ordinary PCs. On a 2.3GHz Linux PC it takes about 30 
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moi 


0.05952 


0.1400 


0.140 


0.150 


mo2 


0.05952 


0.0600 


0.100 


0.050 


b 


1.8412 


1.8400 


10.00 


5.00 


■^ADM 


0.11572 


0.18871 


0.2263 


0.1881 


J Art A f 


0.02315 


0.04233 


0.122 


0.0527 


u> 


0.038 


0.048 


0.0052 


0.013 


di2 


4.224 


4.174 


20.09 


10.20 


di 


2.276 


1.711 


1.70 


1.63 


d2 


2.276 


2.398 


1.96 


2.25 


XCM 





0.74 


1.6 


2.4 


^vncLX 1 


0.0285 


0.106 


0.108 


0.127 


Q.TncLx2 


0.0285 


0.0282 


0.0588 


0.0235 


•^inl 


+0.975 


+1.174 


+9.19 


+4.25 


Xmaxl 


+2.112 


+2.029 


+10.04 


+5.07 


Xoutl 


+3.251 


+2.885 


+10.89 


+5.88 


Xin2 


-0.975 


-0.927 


-9.07 


-4.00 


Xmax2 


-2.112 


-2.145 


-10.05 


-5.13 


Xout2 


-3.251 


-3.325 


-11.03 


-6.25 


Madm (K.o\n/2 

iW(7) ^ Kj ' 


1.7524 


2.8578 


3.427 


2.849 


moi/MADM 


0.5143 


0.7419 


0.619 


0.797 


mQ2/MADM 


0.5143 


0.3179 


0.442 


0.266 


di2/MADM 


36.50 


22.12 


88.78 


54.23 


di/MADM 


19.67 


9.067 


7.51 


8.67 


d2/MADM 


19.67 


12.71 


8.66 


12.0 


loMadm 


0.0044 


0.0091 


0.0012 


0.0024 


Jadm/M^jj]^^ 


1.729 


1.189 


2.382 


1.489 



Table 1. Properties of initial data for different parameters moi, mo2 and b. The 
numbers are first given in units of G = c = k = 1. The total ADM mass Madm is also 
given in solar masses where hq ~ 5 x lO^m^ and n ~ 1 (Note that k^^^^GMadm/c^ is 
dimensionless) . After that we also list some quantities in geometric units {G = c = 1) 
in terms of the total ADM mass. 

hours to push the Hamiltonian constraint down to 10^^ if we use ra^ = = 18 points. 

The main reason for the low memory footprint is that our spectral code needs only very 

few grid points to achieve the quoted accuracies. 

Figure H] shows q in the xy-pla.ne for a binary with rest masses moi = 0.14 and 
= 0.06. We can see how the domain boundaries are adapted such that q is non-zero 

only in the inner domains. Note that q is the rest mass density for k, = n = 1. In 

Figs. [5] and [6] we show the conformal factor ip and the largest shift component By for 

the same configuration. Note that unlike q both are smooth (C^) across the domain 

boundaries. 

In order to further verify our code we have performed a comparison with previous 
results from Taniguchi et al. [21]. In Figs. [7] and Owe show how Madm and Jadm vary 
as a function of uj for a binary with rest masses moi = 0.1461 and mo2 = 0.1299. As 
we can see our results (squares) approach the expected post-Newtonian results (taken 
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Figure 4. q in tlie xy-plane for a binary with rest masses moi = 0.14, mo2 = 0.06 
and b = 1.84. 



Psi 




y 2 




Figure 5. Conformal factor ip in the xy-plane for a binary with rest masses moi = 0.14 
and = 0.06. 



from [371 [251 [38l [39] ) for point particles (dotted line) for small u. At intermediate u our 
results differ from point particle results and instead agree with previous results obtained 
by Taniguchi et al. [21]. Note that, while the agreement in Madm does not look as good 
as for Jadm, the values for Madm from both methods still agree to better than 0.05%. 
With our current code we can construct initial data only up to ~ 0.07. Beyond that 
point, already the first iteration of our elliptic solver fails. We suspect that our initial 
guess of simply using two spherical TOV stars with 5* = is not good enough for close 
configurations, and that the solver would succeed if we provided a guess that is closer to 
the true solution. Notice that while Taniguchi et al. [21] can extend their sequence to 
higher u they also did not record a turning point in either curve for this configuration. 
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Figure 6. Largest shift component By in inertial coordinates for a binary with rest 
masses moi = 0.14 and too2 = 0.06. 
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Figure 7. The ADM mass for a binary with rest masses moi = 0.1461 and 
= 0.1299 as a function of angular velocity in units ofG = c = K = l. Shown 
are results for post-2-Newtonian point particles (dotted line), values from previous 
work [21] (crosses), and results from our new code (squares). 



5. Discussion 

The purpose of this paper is to introduce a new numerical method for the computation 
of binary neutron star initial data with the SGRID code [23]. The method uses six 
domains with different coordinate systems in each domain. The coordinates in four of 
these domains are closely related to the ones suggested in [30|. We have, however, added 
two extra domains with Cartesian coordinates to remove coordinate singularities. All 
our fields are C°° inside each domain. This allows us to use an efficient pseudo-spectral 
collocation method to solve the elliptic equations f|T2l) associated with the initial data 
construction. Note that, we directly solve the 5 Eqs. ( fT2l) . i.e. we do not split the 
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Figure 8. The ADM angular momentum for the same cases as in Fig. [7] in units of 



G = C=K= 1. 



shift in the momentum constraint into a vector and a gradient of a scalar, which would 
introduce an additional elliptic equation. Thus we have one less equation to solve than 
in the original Wilson-Mathews approach |9l |T0] . Since two of our domains extend to 
spatial infinity we are able to easily impose the boundary conditions in Eq. (HM without 
the need of any approximations such as Robin boundary conditions. 

At present we have only considered corotating configurations. The numerical 
method, however, could be easily extended to configurations with arbitrary spins, e.g. by 
following the approach in [20]. This approach involves the addition of one more elliptic 
equation for a velocity potential. The boundary conditions for this extra equation 
need to be imposed at the star surface, which is somewhat involved if one uses cubical 
domains as in [20]. We our new method, however, such boundary conditions can be 
easily imposed, since each star surface is a domain boundary. 

It is well known that the star surfaces of close configurations can develop cusps 
due to tidal forces [191 El EQ]. We have not investigated this issue with our new 
method yet, because our elliptic solver currently fails already during the first step for 
close configurations. We suspect that we need an initial guess that is better than two 
spherical TOV stars with vanishing shift B^. We would like to point out, however, that 
our method should not have any additional problems with such cups if they occur only 
along the x-axis (the line connecting the two stars). The reason is that the a±{B,(j)) 
which appear in the coordinate transformations and describe the star surfaces, can easily 
be chosen such that the domain boundaries have arbitrary cusps on the x-axis. Note that 
such cusp producing a±{B,(j)) are themselves perfectly smooth in A,B,(j) coordinates, 
so that we do not expect to loose spectral accuracy. 
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